﻿using System;
using System.Diagnostics;
using System.Linq;
using System.Threading.Tasks;
using Cudafy;
using Cudafy.Host;
using Cudafy.Translator;

//                         ..                          . .                    
//                         .....................,,,,,,,,,,,.                   
//                         :~::::::::::::::::::::::::::::::,                   
//                      ...::::::::::::::::::::::::::::::::,.                  
//                   ....~?$8888888Z~,,,,,,,,,,::::::::::::,                   
//              ....+,,::IOO=:::~=?IO88Z:,,:::::::::::::::::                   
//             ..~::?D888Z7I::::::?..:+I88I:::::::::::::::::                   
//         ...,~~~Z88O?~:,.:DDD87::::~..~?8D=:::~~~~~~~~~~~:                   
//         .I=+=O88?~..?===OI+?I788O===+..:?DD~============:                   
//      ..I?+?$8D+,.I+??8D88?,...,~?88++?...:I8D++++++???+?:                   
//      .?IIII8$:.III?8DI=:=???O. ..~7I?I8+...=II???I??IIII:.                  
//      .:O7I7OZ..?IIID~...:IIIIO..=IIII8?,..77IIIIIIIII7II:.                  
//      ..~7777O,.:877$$.  :777778I777O8+,..7777788I7777777:.                  
//       ..=7$$$O,.~8$$77..:$$$$$$$$O8I:..7$$$$88I788O$$$$$:.                  
//        ..=8$$$O,.~8$$$$.:$$$$$$D8I~.,I$$$$88+,..:+$8O$$$:                   
//         ..~8$$$77.,+8ZZ$O8888OI=,.I$Z$$$DO+,.  . .,7$$$$:                   
//         . .:7$Z$$O..:?O88=,::,,7$$$ZZZDZ+,.  . .$$ZZZZZ$:                   
//            ..+OZZZZ7=..:=$$ZZZZZZZ88Z?~......Z$ZZZZZZZZZ:                   
//             ..,+$8$ZZ$Z?O888888OI?=:....~$ZZZZZZZZZZZZZZ:.                  
//                ..:?$8DO$O=:::,....+8IZZZZZZZZZZZZZZZZZZZ,.                  
//                  ...,~=?IZZZZZZ$ZZZZZZZZZZZZZZZZZZZZZZZ$,.                  
//                       ..:ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ$,.                  
//                        .:8888888888888888888888888888888,                   
//                        ..,:::::::::::::::::::::::::::::,..                  
//                                                                             

namespace CudafyTsp
{
    internal class Program
    {
        private static void Main()
        {
            new Tsp().FindRoute(Cities.List);
        }
    }

    public class Tsp
    {
        private const int Cities = 10;
        private int _permutations = 1;
        private readonly float[] _latitudes = new float[Cities];
        private readonly float[] _longitudes = new float[Cities];

        public void FindRoute(string cities)
        {
            Console.WriteLine("Starting");
            BuildCities(cities);
            Console.WriteLine("Permutations " + _permutations);
            Console.WriteLine("CpuTsp " + CpuTsp());
            Console.WriteLine("MpuTsp " + MpuTsp());
            Console.WriteLine("GpuTsp " + GpuTsp());
            Console.WriteLine("Done");
            Console.ReadKey(true);
        }

        private void BuildCities(string cities)
        {
            var i = 0;

            foreach (var fields in cities.Split('\r').Select(city => city.Replace("\r", "").Split(',')).Where(fields => fields.Length == 3))
            {
                if (i < Cities)
                {
                    _latitudes[i] = float.Parse(fields[1]);
                    _longitudes[i] = float.Parse(fields[2]);
                    i++;
                }
            }

            for (i = 2; i <= Cities; i++)
            {
                _permutations *= i;
            }
        }

        private Answer CpuTsp()
        {
            var bestDistance = float.MaxValue;
            var bestPermutation = -1;
            var stopWatch = Stopwatch.StartNew();

            for (var permutation = 0; permutation < _permutations; permutation++)
            {
                var path = new int[1, Cities];
                var distance = FindPathDistance(_permutations, permutation, Cities, _latitudes, _longitudes, path, 0);
                if (distance < bestDistance)
                {
                    bestDistance = distance;
                    bestPermutation = permutation;
                }
            }

            return new Answer { Distance = bestDistance, Permutation = bestPermutation, Milliseconds = stopWatch.ElapsedMilliseconds };
        }

        private Answer MpuTsp()
        {
            var bestDistance = float.MaxValue;
            var bestPermutation = -1;
            var locker = new Object();
            var stopWatch = Stopwatch.StartNew();

            Parallel.For(0, _permutations, permutation =>
            {
                var path = new int[1, Cities];
                var distance = FindPathDistance(_permutations, permutation, Cities, _latitudes, _longitudes, path, 0);
                lock (locker)
                {
                    if (distance < bestDistance)
                    {
                        bestDistance = distance;
                        bestPermutation = permutation;
                    }
                }
            });

            return new Answer { Distance = bestDistance, Permutation = bestPermutation, Milliseconds = stopWatch.ElapsedMilliseconds };
        }

        private const int ThreadsPerBlock = 256;
        private const int BlocksPerGrid = 256;

        private Answer GpuTsp()
        {
            gpu.LoadModule(km);

            var stopWatch = Stopwatch.StartNew();
            var answer = new float[BlocksPerGrid, 2];
            var gpuLatitudes = gpu.CopyToDevice(_latitudes.ToArray());
            var gpuLongitudes = gpu.CopyToDevice(_longitudes.ToArray());
            var gpuAnswer = gpu.Allocate(answer);

            gpu.Launch(BlocksPerGrid, ThreadsPerBlock,"GpuFindPathDistance", _permutations, Cities, gpuLatitudes, gpuLongitudes, gpuAnswer);
            gpu.Synchronize();
            gpu.CopyFromDevice(gpuAnswer, answer);

            var bestDistance = float.MaxValue;
            var bestPermutation = 0;
            for (var i = 0; i < BlocksPerGrid; i++)
            {
                if (answer[i, 0] < bestDistance)
                {
                    bestDistance = answer[i, 0];
                    bestPermutation = (int)answer[i, 1];
                }
            }

            return new Answer { Distance = bestDistance, Permutation = bestPermutation, Milliseconds = stopWatch.ElapsedMilliseconds };
        }

        [Cudafy]
        public static void GpuFindPathDistance(GThread thread, int permutations, int cities, float[] latitudes, float[] longitudes, float[,] answer)
        {
            var threadIndex = thread.threadIdx.x; // thread index within the block
            var blockIndex = thread.blockIdx.x; // block index within the grid
            var threadsPerBlock = thread.blockDim.x;
            var blocksPerGrid = thread.gridDim.x;
            var threadsPerGrid = threadsPerBlock * blocksPerGrid;
            var permutation = threadIndex + blockIndex * threadsPerBlock;

            var paths = thread.AllocateShared<int>("path", ThreadsPerBlock, Cities);
            var bestDistances = thread.AllocateShared<float>("dist", ThreadsPerBlock);
            var bestPermutations = thread.AllocateShared<int>("perm", ThreadsPerBlock);

            var bestDistance = float.MaxValue;
            var bestPermutation = 0;

            while (permutation < permutations)
            {
                var distance = FindPathDistance(permutations, permutation, cities, latitudes, longitudes, paths, threadIndex);
                if (distance < bestDistance)
                {
                    bestDistance = distance;
                    bestPermutation = permutation;
                }
                permutation += threadsPerGrid;
            }

            bestDistances[threadIndex] = bestDistance;
            bestPermutations[threadIndex] = bestPermutation;

            thread.SyncThreads();

            for (var i = threadsPerBlock / 2; i > 0; i /= 2) // credit: CUDA By Example, page 79. http://www.amazon.com/CUDA-Example-Introduction-General-Purpose-Programming/dp/0131387685
            {
                if (threadIndex < i)
                {
                    if (bestDistances[threadIndex] > bestDistances[threadIndex + i])
                    {
                        bestDistances[threadIndex] = bestDistances[threadIndex + i];
                        bestPermutations[threadIndex] = bestPermutations[threadIndex + i];
                    }
                }
                thread.SyncThreads();
            }

            if (threadIndex == 0)
            {
                answer[blockIndex, 0] = bestDistances[0];
                answer[blockIndex, 1] = bestPermutations[0];
            }
        }

        [Cudafy]
        public static float FindPathDistance(int permutations, int permutation, int cities, float[] latitudes, float[] longitudes, int[,] paths, int pathIndex)
        {
            PathFromRoutePermutation(permutations, permutation, cities, paths, pathIndex);

            float distance = 0;
            var city = paths[pathIndex, 0];
            var previousLatitude = latitudes[city];
            var previousLongitude = longitudes[city];

            for (var i = 1; i < cities; i++)
            {
                city = paths[pathIndex, i];
                var latitude = latitudes[city];
                var longitude = longitudes[city];
                distance += (float)Math.Sqrt(Square(latitude - previousLatitude) + Square(longitude - previousLongitude));
                previousLatitude = latitude;
                previousLongitude = longitude;
            }

            return distance;
        }

        [Cudafy]
        public static float Square(float value)
        {
            return value * value;
        }

        [Cudafy]
        public static float PathFromRoutePermutation(int permutations, int permutation, int cities, int[,] paths, int pathIndex)
        {
            for (var i = 0; i < cities; i++)
            {
                paths[pathIndex, i] = i;
            }

            for (int remaining = cities, divisor = permutations; remaining > 0; remaining--) // Credit: SpaceRat. http://www.daniweb.com/software-development/cpp/code/274075/all-permutations-non-recursive
            {
                divisor /= remaining;
                var index = (permutation / divisor) % remaining;

                var swap = paths[pathIndex, index];
                paths[pathIndex, index] = paths[pathIndex, remaining - 1];
                paths[pathIndex, remaining - 1] = swap;
            }

            return 0;
        }

        GPGPU gpu = CudafyHost.GetDevice();
        CudafyModule km = CudafyTranslator.Cudafy();
    }

    internal class Answer
    {
        public float Distance;
        public int Permutation;
        public long Milliseconds;

        public override string ToString()
        {
            return Milliseconds + " " + Permutation + " " + Distance;
        }
    }

    internal class Cities
    {
        public static string List =
@"Albany,42.67,73.75
Albuquerque,35.08,106.65
Amarillo,35.18,101.83
Anchorage,61.22,149.9
Atlanta,33.75,84.38
Austin,30.27,97.73
Baker,44.78,117.83
Baltimore,39.3,76.63
Bangor,44.8,68.78
Birmingham,33.5,86.83
Bismarck,46.8,100.78
Boise,43.6,116.22
Boston,42.35,71.08
Buffalo,42.92,78.83
Carlsbad,32.43,104.25
Charleston,32.78,79.93
Charleston,38.35,81.63
Charlotte,35.23,80.83
Cheyenne,41.15,104.87
Chicago,41.83,87.62
Cincinnati,39.13,84.5
Cleveland,41.47,81.62
Columbia,34,81.03
Columbus,40,83.02
Dallas,32.77,96.77
Denver,39.75,105
Des Moines,41.58,93.62
Detroit,42.33,83.05
Dubuque,42.52,90.67
Duluth,46.82,92.08
Eastport,44.9,67
El Centro,32.63,115.55
El Paso,31.77,106.48
Eugene,44.05,123.08
Fargo,46.87,96.8
Flagstaff,35.22,111.68
Fort Worth,32.72,97.32
Fresno,36.73,119.8
Grand Junction,39.08,108.55
Grand Rapids,42.97,85.67
Havre,48.55,109.72
Helena,46.58,112.03
Honolulu,21.3,157.83
Hot Springs,34.52,93.05
Houston,29.75,95.35
Idaho Falls,43.5,112.02
Indianapolis,39.77,86.17
Jackson,32.33,90.2
Jacksonville,30.37,81.67
Juneau,58.3,134.4
Kansas City,39.1,94.58
Key West,24.55,81.8
Knoxville,35.95,83.93
Las Vegas,36.17,115.2
Lewiston,46.4,117.03
Lincoln,40.83,96.67
Long Beach,33.77,118.18
Los Angeles,34.05,118.25
Louisville,38.25,85.77
Manchester,43,71.5
Memphis,35.15,90.05
Miami,25.77,80.2
Milwaukee,43.03,87.92
Minneapolis,44.98,93.23
Mobile,30.7,88.05
Montgomery,32.35,86.3
Montpelier,44.25,72.53
Nashville,36.17,86.78
Newark,40.73,74.17
New Orleans,29.95,90.07
New York,40.78,73.97
Nome,64.42,165.5
Oakland,37.8,122.27
Oklahoma City,35.43,97.47
Omaha,41.25,95.93
Philadelphia,39.95,75.17
Phoenix,33.48,112.07
Pierre,44.37,100.35
Pittsburgh,40.45,79.95
Portland,43.67,70.25
Portland,45.52,122.68
Providence,41.83,71.4
Raleigh,35.77,78.65
Reno,39.5,119.82
Richfield,38.77,112.08
Richmond,37.55,77.48
Roanoke,37.28,79.95
Sacramento,38.58,121.5
Savannah,32.08,81.08
Seattle,47.62,122.33
Shreveport,32.47,93.7
Sitka,57.17,135.25
Spokane,47.67,117.43
Springfield,39.8,89.63
Springfield,42.1,72.57
Springfield,37.22,93.28
Syracuse,43.03,76.13
Tampa,27.95,82.45
Toledo,41.65,83.55
Tulsa,36.15,95.98
Washington,38.88,77.03
Wichita,37.72,97.28
Wilmington,34.23,77.95";
    }
}
